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Abstract - Spatial patterning can be crucially important for understanding the behavior of 
interacting populations. Here we investigate a simple model of parasite and host populations in 
which parasites are random walkers that must come into contact with a host in order to reproduce. 
We focus on the spatial arrangement of parasites around a single host, and we derive using analytics 
and numerical simulations the necessary conditions placed on the parasite fecundity and lifetime 
for the populations long-term survival. We also show that the parasite population can be pushed 
to extinction by a large drift velocity, but, counterintuitively, a small drift velocity generally 
increases the parasite population. 


Introduction. — The dynamics of coupled popula¬ 
tions have been long studied and the topic continues to 
generate much interest from nearly all branches of science 
as a result of the broad applications it offers - from foxes 
and rabbits to forest fires, from chemical kinetics to disease 
control. This interest is sustained partly by innovations 
in analytical and computational tools, and by the recogni¬ 
tion of the crucial roles played by population discreteness 
and various spatial inhomogeneities. Accounting for these 
aspects leads to surprising complexity in the population 
behavior m- Such effects are usually studied within the 
context of the classical predator and prey models intro¬ 
duced by Lotka and Volterra mis], in which the survival 
of each of two competing populations depends directly on 
its interaction with the other. This paradigm, however, 
represents only one of many types of inter-species interac¬ 
tions. Other types of interactions include symbiosis, com¬ 
petition, and coexistence, all of which abound in ecologi¬ 
cal and biological environments. While an ultimate goal of 
ecological models is an understanding of the co-evolution 
of a web of species IMD], our focus here is modest: the 
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parasite-host (PH) type of interaction, inspired by flea in¬ 
festation of household pets. 

The main distinction of PH interactions is that the two 
species do not compete for survival. Instead, parasites 
reproduce only in the presence of a host which provides 
nutrients and breeding ground. Of course, parasites gen¬ 
erally do not kill the host, so as to continue flourishing 
without having to find another host. Apart from fleas 
and pets, other examples of such interactions in nature 
include brood parasitism in birds, blood-sucking parasites 
in mammals and certain types of viruses at the cellular 
level. We also note that intrinsically similar models have 
been proposed to study topics as diverse as self-catalyzing 
reaction-diffusion systems m and economic growth cen¬ 
ters [12]. PH interactions have been specifically investi¬ 
gated in [T3l|T4] , where the system is filled with a uniform 
distribution of hosts before the reactions with the para¬ 
sites take place. By contrast, we focus on a single host, 
stationary or moving uniformly, and study the spatial dis¬ 
tribution of parasites around the host. Systems with mul¬ 
tiple hosts have been considered previously m, but will 
be mentioned only in passing. 

Our model consists of a host and parasites of constant 
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death rate. Parasites must come into contact with the host 
for reproduction. Typically, the life span of the parasites 
is much shorter than that of the host, so we consider only 
the birth-death process of the former, letting our hosts 
be simply immortal. In general, both species are mobile 
and their actions may depend on the location of the oth¬ 
ers. For simplicity, we will begin with a single stationary 
host, while allowing the parasites to perform only random 
walks. In the language of statistical mechanics models, 
this PH model belongs to the class of contact processes 
[16] . in which the spatial structure of each population is 
expected to play significant roles. In our study, the only 
interesting spatial structure is that associated with the 
parasites and we show that this structure is intimately 
linked to the total population, IVtot. In particular, we dis¬ 
cover that, contrary to expectations, A^tot does not mono- 
tonically decrease when the host moves. Formulating this 
system on a discrete lattice, we solve this problem analyt¬ 
ically, with results that agree well with simulations. The 
non-monotonicity of Ntot can be traced to the interplay 
between the biased random walk and the finite carrying 
capacity of the host. 

In the next section, we define the discrete, stochastic 
model in detail, providing a scheme for simulations. We 
then devote the following section to exact theoretical and 
Monte Carlo simulation results for the parasite population 
and distribution. In addition to the analytic solutions, we 
also offer some physical insights, using a continuum ap¬ 
proximation and heuristic arguments. In particular, we 
note that the spatial distribution of parasite takes the 
same form as, say, a “pion cloud” around a nucleon. Fi¬ 
nally we provide a summary and outlook for interesting 
unsolved problems, as well as extensions of this model to 
more realistic behavior of the hosts and the co-evolution 
of the two populations. 

Parasite-host model definition. We consider an 
L‘^ hyper-cubic lattice with periodic boundary conditions. 
Each of its lattice cells, labeled by an integer-valued vec¬ 
tor r, may be occupied by any number of non-interacting 
parasites. The number in a cell at time t is denoted by 
N{r,t). At each time step, each flea dies with probability 
p. The surviving ones then jump to a randomly chosen 
nearest neighbor cell. In our system, there is just a single 
immortal host located at For each parasite in that cell, 
we introduce B new ones and randomly place them in the 
neighboring cellt0. Here, B is the integer part of 

F-V[N{r^,t)] (1) 

where F is the fecundity and V[N{rh,t)], a general Ver- 
hulst factor, models an environment with finite resources 
and depends on the parasite density at r^ ■ For simplicity, 
we will consider only H’s which depend on N/K, where 
K models the host carrying capacity. One consequence is 

^ The offsprings can be placed in the host cell, but that rule leads 
to a large difference between the numbers between the host cell and 
the neighboring ones. 


that K plays the role of setting an overall scale for N, and 
does not enter the general behavior (e.g., extinction) of the 
population. To give F a sensible meaning, we will impose 
H(0) = 1 so that each parasite produces F offsprings in 
the limit oi N K. Although we can analyze the model 
with any H, here, for a variety of reasons m , we use 

V[N] = e-^!^. (2) 

Note that the effective Malthusian growth per time step 
is not just F — p, since death occurs everywhere but birth 
only takes place at t\. 

For sufficiently large E, we expect a steady population 
of parasites: the active state. If there is a single stationary 
host located at = 0, the ensemble average of N{f,t) 
will approach a stationary distribution, p* (r). For a finite 
system with p > 0, there is a finite probability of parasite 
extinction. Since that is an absorbing state, p* = 0 is 
the true stationary state. In this sense, we mean a quasi- 
stationary state when we speak of an active steady state. 

For a non-stationary host, we find more interesting phe¬ 
nomena [15] : however, analytic understanding is challeng¬ 
ing. We focus on a host moving with constant velocity 
along one of the lattice axes (say, x): Fh = vtx. Intu¬ 
itively, one expects the moving host can outrun the para¬ 
sites, so that, for large enough u, the latter goes extinct. 
This picture raises the following questions: Does the par¬ 
asite population always decrease with increasing vl If so, 
what is the critical v for parasite extinction? 

While it is easy to carry out simulations of this model, 
the analysis is less simple. Thus, we turn to a similar 
system: a stationary host with drifting parasites. In this 
case, the parasites perform random walk with a bias e and 
hops to a neighboring cell in direction h, with probability 

l-{h-x)e 
2d ■ 

The statistical properties of this system should be the 
same, under a proper Galilean transformation, as those 
of a host moving with velocity oc ex in a sea of diffusing 
parasites with no drift. Unlike the moving host problem, 
the analysis of drifting parasites is straightforward and, as 
will be presented, simulations of both systems show that 
the differences are indeed minor. 

Without loss of generality, we choose e > 0, i.e., the par¬ 
asites favoring the —x direction, corresponding to the host 
moving along +x with constant velocity v = {v, 0), the sec¬ 
ond entry representing transverse direction. Adopting this 
notation for the rest of this article, we write r = (x,?7), 
and put the host at (0, 0). We implement our agent-based 
simulations on an L x L lattice with periodic boundary 
conditions. For convenience, we choose odd L, so that the 
host is placed at the center. Although we studied various 
lattice sizes, we present mainly the results with L = 101, 
as they provide an adequate picture of the general sys¬ 
tem. With one parasite in each cell initially, we implement 
each Monte Carlo step (MCS) as updating each parasite. 


P-2 





Parasite-host population dynamics 


removing it with probability /r and moving each survivor 
to one of its 2d nearest neighbor cells with probability 
given by eq.([3|). Measuring Nq, the number at the origin, 
we generate B offsprings according to eqs. (fim and place 
each in a randomly chosen nearest neighbor cell. In most 
simulations, we choose K = 100 to provide good statis¬ 
tics and to avoid the absorbing state. We discard the first 
10® MCS to allow the system to reach steady state. We 
then measure N{r) for 4 x 10® MCS after all offsprings 
are placed and before the next culling step. The results 
are used to compile averages at each cell to arrive at the 
stationary distribution 

For the moving host case, the parasites are updated as 
above with e = 0. At the end of the reproductive cycle but 
before the culling steps, we introduce a probability p, with 
which we move the host to the next cell along -\-x. Thus, 
the host velocity is p per MCS. For drifting parasite, the 
average change in x at the end of each MCS is —e/d. We 
compare the two systems with p = e/d. 

We note that the velocities introduced through these 
rules are limited to one lattice spacing per MCS. There 
are many ways to achieve any velocity. However, we do 
not expect further qualitatively novel behavior and will 
restrict ourselves to the rules given above. In the sum¬ 
mary section, we will discuss the physical significance of 
the various parameters. 

Theoretical considerations and simulation re¬ 
sults. — We first study the system with a stationary host 
and parasites drifting along —x. We exploit a mean-field 
equation for the evolution of p{r,t): 

p(r, i + 1) = ^ X! y + a,t)[X + BS{x)S{y + a)] 

a 

+ ^ X! + [-^(1 + + 'r)Hy)] ■ 

T — ±l 

(4) 

Here, A = 1 — /i is the survival probability, <5 is the 
Kronecker delta, and a is one of the 2{d — 1) transverse 
lattice vectors: (0,..., 0, ±1, 0,..., 0). Although B is an in¬ 
teger in simulations, we use it to denote the key quantity 
controlling the birth rate F ■ V = F exp [—po/A'], where 

Po = p(0,0,t). (5) 

The right side of eq.dH) accounts for the occupation at r 
due to transverse hopping of the surviving parasites from 
neighboring cells and possible newborns. In the second 
term on the right side of eq.Q, the effect of the bias is 
incorporated. Except for the non-linear aspect implicit 
in B, solving this equation is straightforward. B only 
depends on p through pQ. Thus, these terms can also be 
written as poBS{...)/2d, so that the only non-linearity is 
contained in a single quantity, po. 

Turning to our main interest, p*{r) in the steady state, 
we defer details of its derivation to the Appendix (which 
also contains the generalization to arbitrary V) and report 


the findings here. Since all of our simulations are on a d = 
2 square lattice, we will drop d and write r = {x,y), etc. 
Using a Fourier transform and regarding the stationary Pq 
as a (to-be-determined) parameter, we find: 

P*ix,y) = PoBG{x,y;X,e) ( 6 ) 


where G is essentially the lattice Green’s function for bi¬ 
ased diffusion with dissipation on a periodic square lattice. 
The unknown can now be fixed: Pq = pQBa{X,e), where 


k,p 


cosfc -I- cosp 

2 — X (cos k -I- cosp) + ieX sin k 


(7) 


embodies the return probability of this particular random 
walker. Here, k and p are integers (1,...,L) times 27r/L. 
a depends not only explicitly on the system size, survival 
and bias, but also implicitly on the boundary conditions. 
With periodic boundary conditions, a is, despite the pres¬ 
ence of e in eq. actually a function of e^. This is hardly 
surprising given that a is analytic in e and cannot depend 
on its sign. As expected, a decreases monotonically as e 
increases for fixed A. In particular, a' = da/d (e^) is neg¬ 
ative at £ = 0. This, however, does not guarantee that the 
total parasite population shares the same behavior. 

Apart from the extinction pg = 0 from pg = PgHcr, the 
non-trivial steady state is given by Ba{X,e) = 1. For the 
special V we chose, the result is 

p*=K\n{Fa) (8) 

while the total parasite population is 

Ntot = — In(Fcr) (9) 

fl<7 

The survival condition is F > F^i^ = l/o'(A,£), rather 
than the naive F > p. The contours in figlTJhelp visualize 
the increasing threshold for Fmin- In the appendix, we 
show that this condition remains valid for any V. 



pai'asite drift, £ 


Fig. 1: Contours of Fmin in the case of L = 101 to sustain a 
population. For example, for F — 4, the parasites goes extinct 
for (A, e) below and to the right of the green contour line. 

We see from eq.® that the bias affects Wot through 
competing factors. As a result, the sign of d^Ntot oc 1 — 
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Fig. 2: (a) F — e phase diagram for A = 0.99, L — 101. The 
green/yellow areas are associated with parasite populations 
larger/smaller than that with no bias. Along the Fa {X,e) = e 
line (dot-dashed), the population is maximal, (b) Ntot{£)/K 
for F — b, corresponding to traversing the solid blue line in 
(a), shows the competing effects of bias. The solid blue line 
is from eq. and the symbols are results from Monte Carlo 
simulations. The agreement is within ~ 1%. 


In(FtT) can change. In particular, for F > e/a{X,0), the 
parasite population first increases as e increases, reach¬ 
ing a maximum when a = e/F, and then decreases to 
extinction at ct = 1/F. For what e will A^tot return to 
the unbiased value? The answer is e{X,F), the solution 
to a (A, 0) ln(FCT (A, e)) = a (A, e) ln(FCT (A, 0)). Since the 
population is enhanced or suppressed on either side of this 
line, it is instructive to show it in the F-e phase diagram. 
We illustrate in figHKa) the case with A = 0.99. We also 
display the line of maximal Wot defined by e, as well as 
the region for extinction. 

In fig|2Kb), we provide an example of the non-monotonic 
behavior corresponding to traversing the F = 5 line (solid 
blue) in fig|2{a). Our simulations confirm the numerical 
predictions. Due to the maximum velocity imposed by 
our rules, the extinction phase exists only for small F. In 
the next section, we show that, for large enough L and 
velocity, extinction prevails for any F as long as /r > 0. 

In addition to Wot, we examine the spatial distribu¬ 
tion of the parasites. As illustrated in figEl a symmetric 
cloud of parasites develops around the host when no bias 
is present. With increasing bias, the cloud is distorted, 
dropping sharply to the right of the host, as new para¬ 
sites drift to the left. These clouds decay exponentially, 
an expected feature associated with death and diffusion, 
with two characteristic length scales as a result of the bias. 
In figlU we show integrated profiles P*(^i v)) of these 
clouds. From these hgures, we see an excellent agreement 
between simulation data and this theory. 

Finally, we turn to the case of a moving host with 
a population of unbiased parasites. The full stochas¬ 
tic process can be written and, naively, there should be 
no difference between this process and the one above af¬ 
ter a Galilean transformation. To compare the two, we 
show the integrated density profiles for a typical case 
{F,X,L = 2.5,0.99,101) with e = 0.5 and p = e/2 in 
figEl The average overall properties of the two system are 
indeed very similar. While the moving host problem may 
be analytically solvable, this result lessens the urgency for 



Fig. 3: Countour plots of simulation data for the density pro¬ 
files of the parasites with drift e = 0,0.5 and 1. In all cases, 
F,X,L = 2.5, 0.99,101 and the host is located at the origin. 


such a pursuit. 

Insights from the continuum limit. To gain 
more intuitive physical insight, we consider the contin¬ 
uum and thermodynamic limits. Here, r is taken to be 
continuous and p{f) is the parasite number density. We 
can write dtp = [DV^ — udx — t“^] p, where D is the dif¬ 
fusion constant, —ux is the drift velocity, and r is the 
average lifetime. Without drift, the characteristic length 
isC = '/Dt. We can cast the effects of the bias in a dimen¬ 
sionless quantity, w = ut j2^. In the absence of a source, 
p vanishes as t —>■ oo. If we place a source of “charge” Q 
(with units of parasite births per unit time) at the origin, 
then a stationary state p*{^ exists and satisfies: 

[_f V" + 2w£,dx + 1] p* (F) = Qr<5(r) (10) 

We see Wot = Qt with an implicit tu-dependence being 
through Q. Notably, for w = 0, p* \s the solution of 
the inhomogeneous Helmholtz equation and, in d = 3, is 
analogous to the potential associated with the Yukawa in¬ 
teraction from particle physics and to the Debye-Hiickel 
interaction from colloidal physics. The behavior far from 
the origin is dominated by e“FI/4 for any d, while it di¬ 
verges as near the origin for d > 2. Unlike in typical 
physics problems, Q here is specihed by the birth rate at 
the host, which is determined by the return probability of 
a newborn to the host. This non-linear feedback produces 
an unusual twist to the physics of finite-ranged interac¬ 
tions, allowing for complex and counterintuitive effects. In 
particular, if there is a true “point” host, then Q depends 
on p(()) and we must deal with ultraviolet singularities d 
for d > 2. The more tractable method is to introduce a 
cutoff, say within radius a, where the reproduction occurs. 

^ Since we consider parasites with r > 0, there are no infrared 
singularities in d < 2. 
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Fig. 4: Integrated density profile of the parasites for both sim¬ 
ulation (symbols) and theory (line) with e = 0, 0.5 and 1. In all 
cases, F,X,L — 2.5, 0.99,101. Parasites of the same horizontal 
coordinates are summed together. 

On the other hand, the precise shape of the reproductive 
region, as long as it is finite, is irrelevant for the large r 
behavior of p. There, the parasite cloud depends only on ^ 
and an “effective charge” Qeff- This approach is similar to 
the renormalization program in quantum and statistical 
field theories, where the large scale properties are univer¬ 
sal, dependent on renormalized quantities (in this case, 
Qeff and and independent of other details of the micro¬ 
scopies. Therefore, the relationship between the details of 
the host and say, iVtot, is not accessible in general. Nev¬ 
ertheless, we illustrate an intuitive and qualitative picture 
using the following arguments. 

We allow reproduction by the parasites that get close 
to the host and denote their numbers by W- Each gives 
birth at a fixed rate $ (fecundity per unit time). We 
also introduce a Verhulst factor, V(JV^), a specific example 
being exp[—Nr/K]. With no bias, we may assume that 
Nr is the integral of p, which is ~ for small r, apart 
from logarithms in d = 2 over a region of size a. We 
can approximate ~ A^tot(a/C)^ and write Q = Nr^V. 
With Ntot = Qt, we arrive at: ~ Nr^V{Nr)T, 

so that a non-trivial solution is 

iVtot~if(a/0“'ln[<i>r(a/ef] . (11) 

Comparing with eq.®, we see that <l’T(a/^)^ plays the 
role of Fa: is the ratio of birth to death rates and 

(a/C)^ is the fraction in the population that can reproduce. 

Now let us look into the effects of biased diffusion. This 
is best revealed in the simple d = 1 case. The parasite 
density profile is 

p{x) = /(^+ -b C-) for x ^ 0 

where (^_) controls the tail along (against) the bias and 
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Fig. 5: Comparison between the integrated density profiles for 
drifting parasites (red o) and moving host (x). Inset gives 
the absolute difference between the two cases with F,X,L = 
2.5,0.99,101. 

+ vFFI. This is similar to the integrated profiles 
shown in figlD Defining W by pdx (a <C ^-), we hnd 
it to be 1 /Vl + times the non-biased case. Though this 
factor decreases monotonically with w, Ntot can increase 
due to the increase in V. The drift reduces the parasite 
population at the host and, in appropriate circumstances, 
allows for more newborns. 

Summary and outlook. — In this article, we report 
studies of a simple system of non-interacting parasites, 
reproducing in the presence of a single host. Using a Ver¬ 
hulst factor for the birth rates, the population reaches 
a steady state at long times, with a non-trivial density 
profile. With a uniformly moving host, the results are es¬ 
sentially the same as those in a system with a stationary 
host and drifting parasites (executing random walks with 
an appropriate bias). This equivalence is fairly intuitive, 
since the two system are related by a Galilean transform 
apart from minor details in stochastic rules. It is tempt¬ 
ing to expect that a moving host is less supportive to the 
parasites, so that their population would decrease with 
the speed. Surprisingly, our study reveals a regime of fe¬ 
cundity in which the total population increases with drift 
speed, before decreasing to extinction. We find analytic 
solutions valid for periodic lattices in any dimension and 
the agreement with simulation is excellent. We also con¬ 
sidered a continuum version to provide a more illuminating 
understanding of the physical system. 

In closing, let us remark on a number of interesting gen¬ 
eralizations of our system, some of which were discovered 
in preliminary studies m- Imposing reflecting EC’s, in¬ 
stead of periodic EC’s, breaks translational invariance and 
leads to novel features: the system supports a larger par¬ 
asite population when the host is placed near a wall or 
a corner due to the increased return probabilities of the 
newborns. An alternative perspective is that the “image 
charge” of the host is closer when it is near a wall or cor¬ 
ner. Similar cooperative behavior emerges when there are 
M > 1 hosts: the total parasite population is more than 
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the sum of the individual cases. As long as the hosts are 
stationary, we can extend our analytic results. 

More interesting properties appear if the host moves, 
either randomly or with deference to the gradient of local 
parasite densities. With reflecting EC’s, it is equally likely 
to find the host anywhere in the former case and more 
likely to find it near the center for the latter. Meanwhile, 
the parasite profile peaks at the corners and the center 
in the two cases, respectively. The details also depend on 
the relative hopping rates, in addition to competing fac¬ 
tors like /r and L. Adding more smart hosts raises the 
natural issue of mutual avoidance. If we probe the tra¬ 
jectories of these hosts, we should find “scattering” events 
because their effective interaction is repulsive, mediated by 
the respective parasite clouds. A systematic exploration, 
which may reveal other novel phenomena, should be un¬ 
dertaken. Here, we limited ourselves to populations with 
relatively low birth rates, so that the systems relax into 
stable steady-states. Regarding our system in the same 
light as the logistic map, Xn+i = Aa;n(l — Xn), we expect, 
for very high birth rates, to encounter instability, possi¬ 
bly transitioning to bifurcation and chaos. Beyond this 
PH system, we can consider multiple species with various 
forms of interdependence, which could pave the way to the 
study of more complex and realistic food-webs. 

Appendix: Solution for the steady state in a fi¬ 
nite discrete lattice. — We present the essentials to de¬ 
rive the stationary solution for eq.(j4]). Writing p*{k,p) = 
'^^^ikx+ip-yp* carrying out the sum on both sides 

of eq.(|4|), we find 


p*{k,p) = pIb 


__ 

d — \A(k, p) + ieX sin k ’ 


( 12 ) 


where A{k,p) = cos + cospi. Inserting this into the 
inverse transform, p*{f) = p where 

p denotes summing over integer (1,..., L) multiples of 
27r/L with a factor of L~‘^, we find eq.® with 


G{x,y;X,e) 


_ ^_ikx-ip-y 

d — XA{k, p) + ieX sin k 

k,p 


(13) 

For self consistency, the unknown pg must satisfy pg = 
p*(0, ()) = PqBG{0, 0; A,£), the d = 2 version being eq.©. 
Finally, the total population is 


iVtot = ^P*(r)=p*(0) 
f 


PqB 

1-A 


* 

£p_ 

pa 


(14) 


It is easy to generalize these results to arbitrary H’s 
that depend on p through (j) = pjK. The non-trivial 
stationary condition then reads H((/'g) = l/{Fa) , while 
Ntot = 4'o{K /= iBK/p)<j)QV{4>Q). Extinction is given 
by Fa = 1/H(0) = 1, while the maximum of A^tot occurs 
at (j)Q, which satisfies i^g = — 91nH|^. and translates to a 

line in the (A,£:) plane through l/a{X,e) = FV{^q). 
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